    // NOTE: All references to clockwise motion are viewed from as
// (x, y) => (0, 1) as indicating going up, (1, 0) as indicating going right
module Rectangulate
    
    let benchmark f loops = 
        let timer = new System.Diagnostics.Stopwatch()
        timer.Start()
        let returnValue =
            for i = 0 to loops do f()
        printfn "Elapsed Time: %i" timer.ElapsedMilliseconds

    let private pointIsInList x list =
        list |> Array.fold (fun acc elem -> acc || (elem = x)) false

    let inline private first  (a, _, _) = a
    let inline private second (_, b, _) = b
    let inline private third  (_, _, c) = c

    let FindLongestAugmentedPath (possibleMatches : int[][]) (vertex : int)= 
        let mutable possibleList = [| [| vertex |] |]
        let mutable iterations = possibleMatches.Length
        let mutable depthRemains = true
        while depthRemains && (iterations > 0) do
            iterations <- iterations - 1
            //let branchesToTest = possibleList |> Array.filter (fun i -> possibleMatches.[i].Length > 0)
                    
            //let permutations =[| for i in 0 .. possibleList.Length - 1 -> for j in 0 .. possibleList.[i].Length - 1 -> (i, j) |]
            let list = possibleList // Save the current instance of this mutable value
            
            // Enumerates all the permutations of matches
            let permutations =
                [| 0 .. list.Length - 1 |] |> Array.map (fun i ->
                    let route = list.[i]
                    if route.Length > 0 then
                        let last = route.[route.Length - 1]

                        // Excludes points already in the route
                        let nextMatches = possibleMatches.[last] |> Array.filter (fun elem -> not (route |> pointIsInList elem))
                     //   printfn "->%d:%A" nextMatches.Length nextMatches

                        if nextMatches.Length > 0 then
                            [| 0 .. nextMatches.Length - 1 |] |> Array.map (fun j -> (i, j))
                        else
                            [| |]
                    else
                        [| |])
                |> Array.filter (fun x -> x.Length > 0)
                |> Array.collect (fun x -> x)

          //  printfn "%A" permutations

            depthRemains <- permutations.Length > 0
            if permutations.Length > 0 then
                // Auguments the routes with next possible matches
                possibleList <- permutations |> Array.map (fun p ->
                    let route = list.[fst p]
                    //printfn "%d;  %A" (fst i) route
                    if route.Length > 0 then
                        let last = route.[route.Length - 1]

                        // Excldes points already in the route
                        let nextMatches = possibleMatches.[last] |> Array.filter (fun elem -> not (route |> pointIsInList elem))
                        Array.append route [| nextMatches.[snd p] |]
                        //Array.append route [| 0 |])
                    else
                        [|  |])
                |> Array.filter (fun x -> x.Length > 0)
        possibleList

    let MaximumMatching (left : int[][]) (right : int[][]) =
        let possibleMatches = Array.append (left |> Array.map (fun i -> i |> Array.map (fun j -> j + left.Length))) right
        let matches =
            [| 0 .. left.Length - 1 |] |> Array.map (fun i -> FindLongestAugmentedPath possibleMatches i |> Array.collect (fun x -> x))
            |> Array.collect (fun x -> x)
        //Brute force
        (*
        let permutations =
            [| 0 .. matches.Length - 1 |] |> Array.map (fun i ->
                [| 0 .. matches.[i].Length|] |> Array.map (fun j -> (i, j)))
            |> Array.collect (fun x -> x)

        
        let ApplyMatching (p : (int*int)[]) (matching : int[][]) =
            fst
            *)
        [| 0 .. left.Length - 1 |] |> Array.map (fun i ->
            [| 0 .. 2 .. matches.Length - 1 |] |> Array.fold (fun acc elem -> // Priortizes the last one
                if (matches.[elem] = i) && (elem < matches.Length - 1) then
                    matches.[elem + 1]
                else
                    acc) -1)
            //|> Array.map (fun i -> i - left.Length)

    let MinimumIndependentSet (left : int[][]) (right : int[][]) =
        let maxMatching = MaximumMatching left right
        [| 0 .. maxMatching.Length - 1 |]
            |> Array.map (fun i -> if maxMatching.[i] < 0 then i else maxMatching.[i])
            |> Array.sort

    let IssueError msg = 
        printfn "%s" msg

    (*
    let bin = System.Random()
    let testMatrixTiles = Array2D.init 4 4 (fun i j -> bin.Next(2)))
    *)
    


    //===========================================================================
    // There are two types of data distictions
    // - Whether it is of tile form or vertex form
    // - Whether it is stored as a matrix or as a list
    //===========================================================================
    // Convertex from tile matrix form to vertex list form
    let MatrixToVertexList (shape : int[,]) =
        [| 0 .. shape.GetLength 0 - 1|] |> Array.map (fun i ->
            [| 0 .. shape.GetLength 1 - 1|] |> Array.filter (fun j -> shape.[i, j] = 1)
            |> Array.map (fun j -> (i, j)))
        //|> Array.filter (fun i -> i.GetLength > 0)
        |> Array.collect (fun x -> x)
            
    // Converts from list form to matrix form
    let ListToMatrix (list : (int*int)[]) =
        let width = list |> Array.map (fun tuple -> fst tuple) |> Array.max
        let height = list |> Array.map (fun tuple -> snd tuple) |> Array.max
        let mutable m = Array2D.zeroCreate (width + 1) (height + 1)
        for i = 0 to list.Length - 1 do
            m.[fst list.[i], snd list.[i]] <- 1
        m
        
    // Converts from tile matrix from to vertex matrix form
    let TileMatrixToVertexMatrix (shape : int[,]) =
        let width = shape.GetLength(0) - 1
        let height = shape.GetLength(1) - 1
        Array2D.init (width + 2) (height + 2) (fun i j ->
            shape.[max (i - 1) 0, max (j - 1) 0] |||
                shape.[min i width, max (j - 1) 0] |||
                shape.[max (i - 1) 0, min j height] |||
                shape.[min i width, min j height]
            ) // 1 if all four corners are 1, 0 otherwise

    // Converts from vertex matrix form back to tile matrix form
    let VertexMatrixToTileMatrix (shape : int[,]) =
        let width = shape.GetLength(0)
        let height = shape.GetLength(1)
        Array2D.init (width - 1) (height - 1) (fun i j ->
            shape.[i, j] * shape.[i + 1, j] * shape.[i, j + 1] * shape.[i + 1, j + 1]
            ) // 1 if all four corners are 1, 0 otherwise
            
    // Prints out a Array2D with x and y coordinates in the expected orientation
    let print2D (ls : int[,]) =
        for y = 0 to (Array2D.length2 ls) - 1 do
            for x = 0 to (Array2D.length1 ls) - 1 do
                printf "%d" ls.[x, y]
            printfn ""

    // Debugging function for printing out a vertex list form
    let printVertexArrayList (list : (int*int)[][]) =
        list |> Array.map (fun i -> printfn "%A\n" <| ListToMatrix i)

    let ListToMatrix2 (list : (int * int * int)[]) =
        let width = list |> Array.map (fun tuple -> first tuple) |> Array.max
        let height = list |> Array.map (fun tuple -> second tuple) |> Array.max
        let mutable m = Array2D.zeroCreate (width + 1) (height + 1)
        for i = 0 to list.Length - 1 do
            m.[first list.[i], second list.[i]] <- third list.[i]
        m


    (*
    let test t =
        let a = GrahamScan <| testShape2
        printf ""
    benchmark (test) 10000
    *)

    let GetInitialPoint vertList fx fy =
        let ystart = vertList |> Array.map (fun t -> fy t) |> Array.min
        let firstrow = vertList |> Array.filter (fun t -> fy t = ystart)
        let xstart = firstrow |> Array.map (fun t -> fx t) |> Array.min
        (firstrow |> Array.filter (fun t -> fx t = xstart)).[0]
            
    // Takes a tile-form matrix <polygon> and returns perimeter as a vertex-form list
    // Returns triplets where the third element is the chaaracterization information
    // Categorizes the vertex based on its location in <polyon>:
    //  - 1: A convex point
    //  - 2: A point on the edge
    //  - 3: A concave point
    //  - 4: Point that lies on the interior of the function
    let TracePerimeter (polygon : int[,]) = 
        let vertexList = polygon |> TileMatrixToVertexMatrix |> MatrixToVertexList
        // rest is the cornerless vertexList, these are the only ones that can overlap for left, right, bottom, top
        let xmax, ymax = polygon.GetLength 0, polygon.GetLength 1
        let rest = vertexList |> Array.filter (fun t -> t <> (0, 0) && t <> (xmax, 0) && t <> (0, ymax) && t <> (xmax, ymax))

        let categorizePoint (point : int * int) =
            polygon.[fst point - 1, snd point - 1] + polygon.[fst point - 1, snd point] +
                polygon.[fst point, snd point - 1] + polygon.[fst point, snd point]

        // So that categorize point doesn't run into index issues
        let categorizeSides test x0 y0 x1 y1 =
            rest |> Array.filter (fun t -> test t)
            |> Array.map (fun t -> fst t, snd t, polygon.[fst t - x0, snd t - y0] + polygon.[fst t - x1, snd t - y1])

        let left = categorizeSides (fun t -> fst t = 0) 0 1 0 0
        let right = categorizeSides (fun t -> fst t = xmax) 1 1 1 0
        let bottom = categorizeSides (fun t -> snd t = 0) 1 0 0 0
        let top = categorizeSides (fun t -> snd t = ymax) 1 1 0 1
        //printfn  "%A\n%A\n%A\n%A" (ListToMatrix2 left) (ListToMatrix2 right) (ListToMatrix2 bottom) (ListToMatrix2 top)

        let interior =
            rest |> Array.filter (fun t -> fst t > 0 && fst t < xmax && snd t > 0 && snd t < ymax)
            |> Array.map (fun t -> fst t, snd t, categorizePoint t)
        
        // Take the only point that left and bottom could possibly share, head, this point must be convex
        let h  = GetInitialPoint vertexList fst snd
        let head = fst h, snd h, 1

        let unorderedPerimeter =
            [| [| (xmax, 0, 1); (0, ymax, 1); (xmax, ymax, 1) |];
                left; right; bottom; top; interior |]
            |> Array.collect (fun x -> x |> Array.filter (fun t ->
                third t <> 4 && (vertexList |> Array.fold (fun acc elem ->
                    acc || (fst elem = first t && snd elem = second t)) false)))
        //let sanitycheck = unorderedPerimeter |> Array.map (fun t ->
        //    (unorderedPerimeter |> Array.filter (fun i -> i = t)).Length) |> Array.max
        //printfn "sanitycheck: %d" sanitycheck
        let head2 = GetInitialPoint unorderedPerimeter first second

        [| 1 .. unorderedPerimeter.Length - 1 |] |> Array.fold (fun acc i ->
            //printfn "%A, %d" acc i
            let last = Array.get acc (i - 1)
            let cur = Array.get acc i
            let lastDir = (first cur - first last, second cur - second last)
            //if third cur = 4 then IssueError "Error: TracePerimeter has entered into the shape"
            let dir = 
                if   third cur = 1 then (-snd lastDir, fst lastDir) // Convex turn counter-clockwise
                elif third cur = 3 then (snd lastDir, -fst lastDir) // Concave, turn clockwise
                else lastDir
            let nextLoc = (first cur + fst dir, second cur + snd dir)
            let next = unorderedPerimeter |> Array.filter (fun t -> first t = fst nextLoc && second t = snd nextLoc)
            //printfn "%A + %A = %A <-> %A:\n%A" cur dir nextLoc next interior
                
            Array.append acc next) [| head; head2 |]
        


    // Requires points <a> to be closer to the origin than <b>,
    // basically that a < b along the axis that they differ
    // Seperate
    let DrawSplitLine (a : int * int * int) (b : int * int * int) (list : (int * int * int)[]) =
        let pA = [| 0 .. list.Length - 1|] |> Array.fold (fun acc elem -> if a = list.[elem] then elem else acc) -1
        let pB = [| 0 .. list.Length - 1|] |> Array.fold (fun acc elem -> if b = list.[elem] then elem else acc) -1
        
        let findOrientation (p1 : int * int * int) (p2 : int * int * int) (p3 : int * int * int) =
           (first p2 - first p1) * (second p3 - second p1) - (second p2 - second p1) * (first p3 - first p1)
        let correctBounds a = (a + list.Length) % list.Length

        let (posA, posB) = if pA < pB then (pA, pB) else (pB, pA)
        //printfn "posA, posB: %A, %A" posA posB
        let chord = // The inbetween points that make up the chord
            if first b - first a = 0 then // Vertical case
                // second a < second b is  always true
                //printfn "Vertical case: %d, %d" (second a) (second b)
                [|second a + 1 .. second b - 1|] |> Array.map (fun i -> (first a, i, 2))
            else // Horizontal Case
                // first a < first b is always true
                //printfn "Horizontal case: %d, %d" (first a) (first b)
                [|first a + 1 .. first b - 1|] |> Array.map (fun i -> (i, second a, 2)) 
        let chordToAdd = if pA < pB then chord else Array.rev chord
        //printfn "chordToAdd: %A" chordToAdd

        // Accounts for the case when pointA is next to pointB (chordToAdd is empty)
        let a = if chordToAdd.Length > 0 then chordToAdd.[0] else list.[posB]
        let b = if chordToAdd.Length > 0 then chordToAdd.[chordToAdd.Length - 1] else list.[posB]

        let aType1 = if findOrientation list.[correctBounds <| posA - 1] list.[posA] a = 0 then 2 else 1
        let bType1 = if findOrientation list.[correctBounds <| posB + 1] list.[posA] b = 0 then 2 else 1
        let aType2 = if findOrientation list.[correctBounds <| posA + 1] list.[posA] a = 0 then 2 else 1
        let bType2 = if findOrientation list.[correctBounds <| posB - 1] list.[posA] b = 0 then 2 else 1
        let point i o = [|(first list.[i], second list.[i], o)|]
            
        let rect1 =
            [| list.[0 .. posA - 1]; point posA aType1; chordToAdd; point posB bType1; list.[posB + 1 .. list.Length - 1] |]
            |> Array.collect (fun x -> (*printfn "%A" x;*) x)
        //printfn ""
        let rect2 =
            [|point posA aType2; list.[posA + 1 .. posB - 1]; point posB bType2; Array.rev chordToAdd|]
            |> Array.collect (fun x -> (*printfn "%A" x;*) x)
        [| rect1; rect2 |]
            

    // Defines a function that will check if a point is within a list
    let private pointIsInList2 (list : (int*int*int)[]) (point : (int*int*int)) =
        list |> Array.fold (fun acc p -> acc || (first p = first point && second p = second point)) false

    let private pointIsInRectangle (rectangles : (int*int*int)[][]) (point : int*int*int) =
        let matches = rectangles |> Array.map (fun i -> pointIsInList2 i point)
        [| 0 .. rectangles.Length - 1|] |> Array.filter (fun i -> matches.[i])





    
    let ChordlessRectangulate (perimeter : (int * int * int)[]) = 
        let concaveVertexIndices = [| 0 .. perimeter.Length - 1|] |> Array.filter (fun i -> third perimeter.[i] = 3)
        let concaveVertices = concaveVertexIndices |> Array.map (fun i -> perimeter.[i])
        
        // Pairs each concave vertex with the direction along which the rectangle will be split
        // the first point should never be a concave vertex, so you can always do index - 1 within limits
        let concaveVertexMatches = concaveVertexIndices |> Array.map (fun i ->
            (first perimeter.[i] - first perimeter.[i - 1], second perimeter.[i] - second perimeter.[i - 1]))
        //printfn "%A" concaveVertexMatches
        
        let mutable rectList = [| perimeter |]
        for i = 0 to concaveVertices.Length - 1 do
            let p = (pointIsInRectangle rectList concaveVertices.[i]).[0]
            let start = concaveVertices.[i]
            let dir = concaveVertexMatches.[i]
            let till = 
                if fst dir = 0 then // Vertical case
                    rectList.[p]
                    |> Array.filter (fun t -> first t = first start)
                    |> Array.filter (fun t -> if snd dir > 0 then second t - second start > 0 else second t - second start < 0)
                    |> Array.sortBy (fun t -> abs <| second t - second start)
                else // Horizontal case
                    rectList.[p]
                    |> Array.filter (fun t -> second t = second start)
                    |> Array.filter (fun t -> if fst dir > 0 then first t - first start > 0 else first t - first start < 0)
                    |> Array.sortBy (fun t -> abs <| first t - first start)
            // Sort along the differing axis for use by DrawSplitLine
            let (a, b) = if fst dir >= 0 && snd dir >= 0 then (start, till.[0]) else (till.[0], start)
            // And replace the pth entry with the two rects from splitRect
            rectList <- Array.append rectList (DrawSplitLine a b rectList.[p])
            rectList <- Array.append rectList.[0 .. p - 1] rectList.[p + 1 .. rectList.Length - 1]
            //printfn "%A" rectList
        rectList |> Array.map (fun r -> r |> Array.map (fun t -> first t, second t))
    
    
    
    // Splits the tile-form shape <tiledShape> into chordless rectangles
    // Chords are when two concave vertices are
    // Note: Error when two chords intersect but the maximum independent set works around this
    // TODO: test everything as done as a sequence
    let Rectangulate (tiledShape : int[,]) =
        let rect : (int*int*int)[] = TracePerimeter tiledShape
        let concaveVertices = rect |> Array.filter (fun t -> third t = 3)
        //printfn "ConcaveVertices: %A" concaveVertices
        
        // Searches for chords along a certain axis orientated by the values give
        // f0 and v0 are the axis shared, f1 and v1 are the axis by which they differ
        let applySearch (list : (int * int * int)[]) f0 f1 v0 v1 = 
            // Some helper functions for triplets
            let inline combine m0 v0 m1 v1 = m0 * first v0 + m1 * first v1, m0 * second v0 + m1 * second v1, 0
            
            [| 0 .. list.Length - 1 |]  |> Array.map (fun f -> // Advance through all the points with triangle pairs
                let from = list.[f] // Note: list is all concave points
                let t =
                    list.[f + 1 .. list.Length - 1]               // For every point after f,
                    |> Array.filter (fun till -> f0 from = f0 till) // find points along the same axis,
                    |> Array.filter (fun till ->                    // remove points with parts of the perimeter in between
                        let axis = f0 from // uses f0 hence same for <from> and <till>
                        let (a, b) = if f1 from < f1 till then f1 from, f1 till else f1 till, f1 from
                        [| a + 1 .. b - 1 |] |> Array.fold (fun acc i ->
                            acc && not (pointIsInList2 rect <| combine axis v0 i v1)) true)
                    
                // Pair the two points together and sort along the axis differ by for DrawSplitLine use
                Array.append [| from |] t |> Array.sortBy (fun tuple -> f1 tuple)) 
            |> Array.filter (fun x -> x.Length = 2) // If there was no match, then it's left over as a empty list
        let horizontal = applySearch concaveVertices second first (0, 1, 0) (1, 0, 0)
        let vertical = applySearch concaveVertices first second (1, 0, 0) (0, 1, 0)
//        printfn "horizontal: %A\nvertical: %A" horizontal vertical
        
        // Finds which chords in list2 intersect with the chords in list
        let inline findMatches (list : ((int * int * int))[][]) (list2 : ((int * int * int))[][]) f0 f1 =
            list |> Array.map (fun tupleList ->
                let w, a, b = f0 tupleList.[0], f1 tupleList.[0], f1 tupleList.[1]
                
                [| 0 .. list2.Length - 1 |] |> Array.filter (fun i ->
                    // Since the chord points are sorted for DrawSplitLine, a < b and c < d
                    let u, c, d = f1 list2.[i].[0], f0 list2.[i].[0], f0 list2.[i].[1]
                    ((a < u) && (u < b) && (c < w) && (w < d))))
        let matchesH = findMatches horizontal vertical second first
        let matchesV = findMatches vertical horizontal first second
        //printfn "matchesH: %A\nmatchesV: %A" matchesH matchesV
            
        // TODO: check if this works in the case of no intersecting chords
        let independentSet = MinimumIndependentSet matchesH matchesV
        let chordsToRectangulateOver = independentSet |> Array.map (fun i -> (Array.append horizontal vertical).[i])
        //printfn "indepedentSet: %A" independentSet
        //printfn "chordsToRectangulateOver: %A" chordsToRectangulateOver
        
        let mutable rectList : (int * int * int)[][] = [| rect |]
        for i = 0 to chordsToRectangulateOver.Length - 1 do
            let chord = chordsToRectangulateOver.[i]
            let a, b = chord.[0], chord.[1]

            let rList = rectList // Get a immutable snapshot of the current rectList
            // Find the index of the rectangle that contains both a and b
            let p = (-1, pointIsInRectangle rList a) ||> Array.fold (fun acc elem ->
                if (pointIsInRectangle rList b |> Array.fold (fun acc2 elem2 -> acc2 || (elem2 = elem)) false)
                then elem else acc)
            rectList <- Array.append rectList (DrawSplitLine a b rectList.[p])
            rectList <- Array.append rectList.[0 .. p - 1] rectList.[p + 1 .. rectList.Length - 1]

        rectList |> Array.map (fun r -> r |> ChordlessRectangulate ) |> Array.collect (fun x -> x)

        
    // Checks orientation of perimeter and if everything is spaced one apart from the previous and next
    let checkPerimeter (list : (int * int)[]) =
        let correctBounds a = (a + list.Length) % list.Length
        let getDir a b = abs (fst a - fst b), abs (snd a - snd b)
        let ystart = list |> Array.map (fun t -> snd t) |> Array.min
        let hlist = [| 0 .. list.Length - 1 |] |> Array.filter (fun i -> snd list.[i] = ystart)
        let xstart = hlist |> Array.map (fun t -> fst list.[t]) |> Array.min
        let head = (hlist |> Array.filter (fun i -> fst list.[i] = xstart)).[0]
        
        let integrity =
            [| 0 ..  list.Length - 1 |] |> Array.fold (fun acc i ->
            let dir1 = getDir list.[correctBounds (i-1)] list.[i]
            let dir2 = getDir list.[correctBounds (i+1)] list.[i]
            acc && fst dir1 + snd dir1 = 1 && fst dir2 + snd dir2 = 1) true
        let orientation = fst list.[correctBounds (head + 1)] - fst list.[head] = 1
        if integrity = false then printfn "Elements are not spaced next to each other"
        if orientation = false then printfn "Some elements are ordered incorrectly"
        integrity && orientation
        (*match perimeter with
        | :? ((int * int * int)[]) as list -> printfn "ohai"
        | :? ((int * int)[]) as list -> printfn "does stuff"
        | _ -> IssueError "Incorrect input to checkPerimeter"; false*)

    let checkPerimeter2 (list : (int * int * int)[]) =
        list |> Array.map (fun t -> first t, second t) |> checkPerimeter

        
    // Debugging function

    let checkRectangulate (rectList : (int * int)[][]) =
        // Test 1 - Each recatangle ordered clockwise and spaced one apart
        if rectList |> Array.exists (fun r -> not <| checkPerimeter r) then false else
        // Test 2 - There's at least 4 points, (one tile means four vertices)
        if rectList |> Array.exists (fun r -> r.Length < 4) then IssueError "Less than 4 vertices"; false else
        
        // Test 3 - Checks for any concave points which rectangle cannot have
        let findOrientation (p1 : int * int) (p2 : int * int) (p3 : int * int) =
           (fst p2 - fst p1) * (snd p3 - snd p1) - (snd p2 - snd p1) * (fst p3 - fst p1)
        
        // Categorize each point wrapping indexed i based on the orientation formed from i-1, i, and i+1
        let catRectList = rectList |> Array.map (fun r ->
            let initialdir = fst r.[1] - fst r.[0], snd r.[1] - snd r.[0]
            let p = [| [| r.[r.Length - 1] |]; r; [| r.[0] |] |] |> Array.collect (fun x -> x)
            [| 1 .. r.Length |] |> Array.map (fun i ->
                let o = findOrientation p.[i - 1] p.[i] p.[i + 1]
                fst r.[i - 1], snd r.[i - 1], if o > 0 then 1 elif o = 0 then 2 else 3))
        if catRectList |> Array.exists (fun r -> r |> Array.exists (fun (_, _, c) -> c = 3))
            then IssueError "List contains a non-rectangle"; false else

        true

    let reconstructRectangulation (perimeterList : (int * int)[][]) = 
        // Gets the bottom-left corner and top-right corner of each rect
        let dimensions = perimeterList |> Array.map (fun r ->
            let x = r |> Array.sortBy (fun t -> fst t)
            let y = r |> Array.sortBy (fun t -> snd t)
            [[fst x.[0]; snd y.[0]]; [fst x.[x.Length - 1]; snd y.[y.Length - 1]]])
        
        // The max of those yield the total size required for the entire rectangle
        let width =  dimensions |> Array.map (fun x -> x.[1].[0]) |> Array.max
        let height = dimensions |> Array.map (fun x -> x.[1].[1]) |> Array.max
        
        // Produces a filled rectangle from the bottom
        // additionally converts from vertex form to tile form
        let rectList =
            [| 0 .. perimeterList.Length - 1 |]
            |> Array.map (fun i ->
                let min, max = dimensions.[i].[0], dimensions.[i].[1]
                [| min.[0] .. max.[0] - 1|]
                    |> Array.map (fun x -> [| for y = min.[1] to max.[1] - 1 do yield x, y |])
                    |> Array.collect (fun x -> x))
        let m = Array2D.zeroCreate width height
        for i = 0 to rectList.Length - 1 do
            //printfn "%A\n" <| ListToMatrix perimeterList.[i]
            for j = 0 to rectList.[i].Length - 1 do
                m.[fst rectList.[i].[j], snd rectList.[i].[j]] <- i + 1
        m

    let ConstructHolelessRectangle (width : int) (height : int) =
        printfn "Place holder"
    (*
    let testShape2 = [| (0,2); (0,3); (1,1); (1,2); (2,1); (2,2); (3,1); (3,2); (4,2);
                        (1, 3); (2,3); (3,3); (4,3); (4,4); (4,5); (5,2);
                        (5,3); (5,4); (5,5); (5,6);
                        (0,6); (0,7); (6,3); (6,4); (6,5);
                        (0,8); (1,8); (2,8); (3,8); (4,8); (7,3); (7,4);
                        (5,8); (5,7);
                        (6,7); (1,6); (1,7); (2,6); (2,7); (3,6); (3,7); (4,6); (4,7);
                        (6,6); (7,6); (7,5); (8,5); (8,4); (8,3); (8,2); (7,2);
                        (6,2); 
                        (6,1); (5,1); (4,1);
                        (4,0); (3,0); (2,0); (1,0); (0,0); (0,1) |]
    printfn "%A" <| ListToMatrix testShape2
    printfn "%A" <| ListToMatrix2 (TracePerimeter <| ListToMatrix testShape2)
    DrawSplitLine (4, 6, 3) (8, 6, 3) (TracePerimeter <| ListToMatrix testShape2)
    
    let chk = Rectangulate <| ListToMatrix testShape2
    chk |> Array.map (fun r -> printfn "%A\n"  <| ListToMatrix r)
    chk |> Array.map (fun r -> r |> Array.map (fun t -> fst t, snd t, 0) |> checkPerimeter)
    checkRectangulate chk |> Array.map (fun r -> ListToMatrix2)

    Rectangulate <| ListToMatrix testShape2
    (Rectangulate <| ListToMatrix testShape2) |> Array.map (fun r -> printfn "%A" <| ListToMatrix r);;
    *)
    (*
    let test t =
        let a = Rectangulate <| ListToMatrix testShape2
        printf ""
    benchmark (test) 20000
    *)